\documentclass{article}
\usepackage{geometry}
\geometry{margin=1.5cm, vmargin={0pt,1cm}}
\setlength{\topmargin}{-1cm}
\setlength{\paperheight}{29.7cm}
\setlength{\textheight}{25.3cm}

% useful packages.
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{enumerate}
\usepackage{graphicx}
\usepackage{multicol}
\usepackage{fancyhdr}
\usepackage{layout}
%\usepackage{ctex}
\usepackage{listings}
\usepackage{subfigure}


% some common command
\newcommand{\dif}{\mathrm{d}}
\newcommand{\avg}[1]{\left\langle #1 \right\rangle}
\newcommand{\difFrac}[2]{\frac{\dif #1}{\dif #2}}
\newcommand{\pdfFrac}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\OFL}{\mathrm{OFL}}
\newcommand{\UFL}{\mathrm{UFL}}
\newcommand{\fl}{\mathrm{fl}}
\newcommand{\op}{\odot}
\newcommand{\Eabs}{E_{\mathrm{abs}}}
\newcommand{\Erel}{E_{\mathrm{rel}}}

\usepackage{xcolor}
\usepackage{fontspec} 
\definecolor{dkgreen}{rgb}{0,0.6,0}
\definecolor{gray}{rgb}{0.5,0.5,0.5}
\definecolor{comment}{rgb}{0.56,0.64,0.68}

\newfontfamily\monaco{Monaco}
\lstset {
aboveskip=3mm,
belowskip=3mm,
showstringspaces=false,       % underline spaces within strings
columns=flexible,
framerule=1pt,
rulecolor=\color{gray!35},
backgroundcolor=\color{gray!5},
basicstyle={\small\monaco},           % the size of the fonts that are used for the code
numbers=left,                   % where to put the line-numbers
numberstyle=\tiny\monaco\color{gray},  % the style that is used for the line-numbers
numbersep=5pt,                  % how far the line-numbers are from the code
commentstyle=\color{comment},
keywordstyle=\color{blue},
stringstyle=\color{dkgreen},
tabsize=2,                      % sets default tabsize to 2 spaces
captionpos=b,                   % sets the caption-position to bottom
breaklines=true,                % sets automatic line breaking
breakatwhitespace=false,        % sets if automatic breaks should only happen at whitespace
escapeinside={\%*}{*)},            % if add LaTeX within your code
morekeywords={*,...}               % if add more keywords to the set
}


\begin{document}
\title{Homework \#6}
\pagestyle{fancy}
\lhead{Name Li HuiTeng 3180102114}
\chead{ numPDE\#6}
\rhead{Date 21.06.20}


\tableofcontents

\newpage

\section{12.5 Generate the fifth order interpolation matrix and explain the coincidence.}
\begin{proof}[Proof]
	The procedures make no difference to 4th order case. Denote $f':=\difFrac{f}{x_1}$ and the only changes is just the following calculation.
	\begin{align*}
		\begin{array}{l}\begin{bmatrix}1&1&1&1&1\\-1&1&-1&1&-1\\2&4&8&16&32\\-2&4&-8&16&-32\\3&9&27&81&243\end{bmatrix}\begin{bmatrix}h&&&&\\&\textstyle\frac{h^2}2&&&\\&&\textstyle\frac{h^3}6&&\\&&&\textstyle\frac{h^4}{24}&\\&&&&\textstyle\frac{h^5}{120}\end{bmatrix}{\begin{bmatrix}\Phi'\\\Phi''\\\Phi'''\\\Phi''''\\\Phi'''''\end{bmatrix}}_{x_1=\overline x}=\begin{bmatrix}\delta_{i+e^1}\\-\delta_i\\\delta_{i+e^1}+\delta_{i+2e^1}\\-\delta_{i-e^1}-\delta_i\\\delta_{i+3e^1}+\delta_{i+2e^1}+\delta_{i+e^1}\end{bmatrix}+O(h^6)\Rightarrow\\{\begin{bmatrix}h\Phi'\\h^2\Phi''\\h^3\Phi'''\\h^4\Phi''''\\h^5\Phi'''''\end{bmatrix}}_{x_1=\overline x}=C\begin{bmatrix}\delta_{i+e^1}\\-\delta_i\\\delta_{i+e^1}+\delta_{i+2e^1}\\-\delta_{i-e^1}-\delta_i\\\delta_{i+3e^1}+\delta_{i+2e^1}+\delta_{i+e^1}\end{bmatrix}+O(h^6).\\\\M=\begin{bmatrix}1&&&&\\&-1&&&\\-1&&1&&\\&1&&-1&\\&&-1&&1\end{bmatrix}\Rightarrow{\begin{bmatrix}\phi\\h\phi'\\h^2\phi''\\h^3\phi'''\\h^4\phi''''\end{bmatrix}}_{x_1=\overline x}={\textstyle\frac1h}T^{(5)}\begin{bmatrix}\delta_{i+e^1}\\\delta_i\\\delta_{i+2e^1}\\\delta_{i-e^1}\\\delta_{i+3e^1}\end{bmatrix}+O(h^5).\\\\\end{array}
	\end{align*}
	Then we compute $C$ and $T^{(5)}$ by Matlab R2019b.
	\lstset{language=Matlab}
	\begin{lstlisting}
%This is a file named pro_125.m.
A=[1 1 1 1 1;-1 1 -1 1 -1;2 4 8 16 32;-2 4 -8 16 -32;3 9 27 81 243];
for i=2:5
A(:,i)=A(:,i)./factorial(i);
end
format rat
C=inv(A)
M=[1 0 0 0 0;0 -1 0 0 0;-1 0 1 0 0;0 1 0 -1 0;0 0 -1 0 1];
format rat
T5=C*inv(M)

%This is the output.
>> pro_125

C =

       1             -1/2           -1/4            1/20           1/30    
       4/3            4/3           -1/12          -1/12           0       
      -7/2           -1/4            7/4           -1/4           -1/4     
      -4             -4              1              1              0       
      10              5             -5             -1              1       


T5 =

      47/60           9/20         -13/60          -1/20           1/30    
       5/4           -5/4           -1/12           1/12           0       
      -2              1/2            3/2            1/4           -1/4     
      -3              3              1             -1              0       
       6             -4             -4              1              1  
\end{lstlisting}
	We can notice that the second row of $T^{(5)}$ is as same as $T^{(4)}$'s.
	To begin with, our $<\phi'>$ formula of $T^{(4)}$ seemingly should have been third-order, since
	\begin{align*}
		h\phi'(\overline{x})=\frac{1}{h}(T^{(4)}\vec{\delta})_2+O(h^4) \Rightarrow \phi'(\overline{x})&=\frac{1}{h^2}(T^{(4)}\vec{\delta})_2+O(h^3) \Rightarrow	\\
		<\difFrac{\phi}{x_d}>_{i+\frac{1}{2}e^d}&=\frac{1}{h}(\frac{5}{4}<\phi>_{i+e^d}-\frac{5}{4}<\phi>_i-\frac{1}{12}<\phi>_{i+2e^d}+\frac{1}{12}<\phi>_{i-e^d}))+O(h^3).
	\end{align*}
	However, the remainder is indeed at least fourth-order. 
	This is because the cancellation from the symmetry of the above difference 
	stencil where all odd order terms cancel each other out and fourth order terms remain. 

	This fact might also be contradictory to our instinct. Since the relation between terms with same coefficient is 'minus' rather than 'add', 
	we were taught in Numerical Analysis that all even order terms should cancel out each other and odd order term should remain, which is just 
	the opposite. 
	To explain it, let's take 1D stencil as an example. Consider $E(x)=\int_\xi^{x}\phi$ where $\xi$ is fixed to $x_{i+1}$ to make $E(x_{i+1})=0$.
	Then the sentence can be reduced to 
	\begin{align*}
		E''(x_{i+1})&=\frac{1}{h^2}(\frac{5}{4}E(x_{i+2})-\frac{-5}{4}E(x_{i})-\frac{1}{12}(E(x_{i+3})-E(x_{i+2}))+\frac{1}{12}(E(x_i)-E(x_{i-1})))+O(h^3)\\
					&=\frac{1}{h^2}(-\frac{1}{12}E(x_{i+3})-\frac{1}{12}E(x_{i-1})+\frac{4}{3}E(x_{i+2})+\frac{4}{3}E(x_{i}))+O(h^3).	
	\end{align*} 
	Under the suppose of $E(x_{i+1})=0$, one can quickly figure out that it's indeed a central difference formula of order 4, with all odd order terms vanishing.

	Hence, the symmetry of the previous difference stencil for $<\phi>$ ensures the absence of all odd order terms and the accuracy of order 4.
	
	Similarly, for $T^{(5)}$,
	the asymmetry forces its maxium accuracy for $\phi'$ to just stay fourth-order. These facts, combining with the uniqueness of interpolation, 
	ensure such a coincidence.
\end{proof}

\section{12.10 Derive the fifth order formulas for ghost cell.}
\begin{proof}[Proof]
	First we have
	\begin{align*}
		 & -<\psi>_{i+e^d}+5<\psi>_{i}-10<\psi>_{i-e^d}+10<\psi>_{i-2e^d}-5<\psi>_{i-3e^d}+<\psi>_{i-4e^d}                                    \\
		 & =\frac{1}{h^D}\int_{C_i}-\psi(x+e^d)+5\psi(x)-10\psi(x-e^d)+10\psi(x-2e^d)-5\psi(x-3e^d)+\psi(x-4e^d):=\frac{1}{h^D}\int_{C_i}S_1.
	\end{align*}
	Denote $f':=\difFrac{f}{x_d}$. Applying Taylor expansions at $x$ repeatedly leads to
	\begin{align*}
		S_1 & =5\phi-(\phi+h\phi'+\frac{h^2}{2}\phi''+\frac{h^3}{6}\phi'''+\frac{h^4}{24}\phi'''')          \\
		    & -10(\phi-h\phi'+\frac{h^2}{2}\phi''-\frac{h^3}{6}\phi'''+\frac{h^4}{24}\phi'''')              \\
		    & +10(\phi-2h\phi'+4\frac{h^2}{2}\phi''-8\frac{h^3}{6}\phi'''+16\frac{h^4}{24}\phi'''')         \\
		    & -5(\phi-3h\phi'+9\frac{h^2}{2}\phi''-27\frac{h^3}{6}\phi'''+81\frac{h^4}{24}\phi'''')         \\
		    & +(\phi-4h\phi'+16\frac{h^2}{2}\phi''-64\frac{h^3}{6}\phi'''+128\frac{h^4}{24}\phi'''')+O(h^5) \\
		    & =O(h^5).
	\end{align*}
	Thus,
	\begin{align*}
		\frac{1}{h^D}\int_{C_i}S_1=\frac{1}{h^D}\int_{C_i}O(h^5)=O(h^5),
	\end{align*}
	which completes the proof.

	For (12.32), denote 
	\[ 
		Sum_2    :=-<\psi>_{i+\frac{1}{2}e^d}+\frac{1}{60}(137<\psi>_{i}-163<\psi>_{i-e^d}+137<\psi>_{i-2e^d}-63<\psi>_{i-3e^d}+12<\psi>_{i-4e^d}). 	
	\]
	Combining with the fifth-order formula derived in Exercise 12.5 and a shift of the above conclusion, we have
	\begin{align*}
		Sum_2   & =\frac{1}{60}(-47<\phi>_{i+e^d}+13<\phi>_{i+2e^d}-2<\phi>_{i+3e^d}+110<\psi>_{i}-160<\psi>_{i-e^d}+137<\psi>_{i-2e^d}             \\
		      & -63<\psi>_{i-3e^d}+12<\psi>_{i-4e^d}+O(h^5)),                                                                                     \\
		60Sum_2 & =12(-<\psi>_{i+e^d}+5<\psi>_{i}-10<\psi>_{i-e^d}+10<\psi>_{i-2e^d}-5<\psi>_{i-3e^d}+<\psi>_{i-4e^d})                              \\
		      & -3(-<\psi>_{i+2e^d}+5<\psi>_{i+e^d}-10<\psi>_{i}+10<\psi>_{i-e^d}-5<\psi>_{i-2e^d}+<\psi>_{i-3e^d})                               \\
		      & +2(-<\psi>_{i+3e^d}+5<\psi>_{i+2e^d}-10<\psi>_{i+e^d}+10<\psi>_{i}-5<\psi>_{i-e^d}+<\psi>_{i-2e^d})+O(h^5)                        \\
		      & =(12-3+2)O(h^5)+O(h^5)=O(h^5),
	\end{align*}
	which completes the proof.
\end{proof}

\section{12.20 Prove Lemma 12.19.}
\begin{proof}[Proof]
	From (12.33), constant u and <f>=0 ensure
	\begin{align*}
		\difFrac{<\phi>_i}{t}=-\frac{1}{h}\sum_{d=1}^Du_d(<\phi>_{i+\frac{1}{2}e^d}-<\phi>_{i-\frac{1}{2}e^d})+\nu L<\phi>_i.
	\end{align*}
	For advection terms, we have
	\begin{align*}
		 & \quad-\frac{1}{h}\sum_{d=1}^Du_d(<\phi>_{i+\frac{1}{2}e^d}-<\phi>_{i-\frac{1}{2}e^d})                 \\
		 & =-\frac{1}{12h}\sum_{d=1}^Du_d(-<\phi>_{i+2e^d}+8<\phi>_{i+e^d}-8<\phi>_{i-e^d}+<\phi>_{i-2e^d})      \\
		 & =-\frac{<\phi>_i}{12h}\sum_{d=1}^Du_d(-e^{2\theta_di}+8e^{\theta_di}-8e^{-\theta_di}-e^{-2\theta_di}) \\
		 & =-\frac{<\phi>_i}{12h}\sum_{d=1}^D2i\cdot u_d(8sin\theta_d-sin2\theta_d)                              \\
		 & =-\frac{<\phi>_i}{12h}\sum_{d=1}^D4i\cdot u_dsin\theta_d(3+2sin^2\frac{\theta_d}{2})                  \\
		 & =-i\frac{<\phi>_i}{h}\sum_{d=1}^Du_dsin\theta_d(1+\frac{2}{3}sin^2\frac{\theta_d}{2})                 \\
		 & :=i\lambda^a<\phi>_i.
	\end{align*}
	For diffusion term, we have
	\begin{align*}
		 & \quad \nu L<\phi>_i                                                                                         \\
		 & =\frac{1}{12h^2}\sum_{d=1}^D-<\phi>_{i+2e^d}+16<\phi>_{i+e^d}-30<\phi>_{i}+16<\phi>_{i-e^d}-<\phi>_{i-2e^d} \\
		 & =\frac{<\phi>_i}{12h^2}\sum_{d=1}^D(-e^{2\theta_di}+16e^{\theta_di}+16e^{-\theta_di}-e^{-2\theta_di}-30)    \\
		 & =\frac{<\phi>_i}{12h^2}\sum_{d=1}^D(-2(1-2sin^2\theta_d)+32(1-2sin^2\frac{\theta_d}{2})-30)                 \\
		 & =-\frac{<\phi>_i}{12h^2}\sum_{d=1}^D(16sin^2\frac{\theta_d}{2}(4+cos^2\frac{\theta_d}{2}))                  \\
		 & =-4\frac{<\phi>_i}{h^2}\sum_{d=1}^Dsin^2\frac{\theta_d}{2}(1+\frac{1}{3}sin^2\frac{\theta_d}{2})            \\
		 & :=\lambda^d<\phi>_i.
	\end{align*}
	Thus, combining the above results leads to
	\begin{align*}
		\difFrac{<\phi>_i}{t} & =(\lambda^d+i\lambda^a)<\phi>_i,           \\
		\difFrac{y}{t}        & =(\lambda^d+i\lambda^a)y(t):=\lambda y(t),
	\end{align*}
	which completes the proof.
\end{proof}

\section{12.21 Reproduce the stability region and determine the bound of stable Courant number.}
\begin{proof}[Proof]
	The plots are reproduced by MATLAB R2019b.
	\lstset{language=Matlab}
	\begin{lstlisting}
%This is a file named ARK4_coeffi.m.
b=zeros(1,6);
aE=zeros(6,6);
aI=zeros(6,6);
gama=0.25;

b(1) = 0.15791629516167136;
b(2) = 0.;
b(3) = 0.18675894052400077;
b(4) = 0.6805652953093346;
b(5) = -0.27524053099500667;
b(6)=gama;

aE(3,1) =0.221776;
aE(3,2) = 0.110224;
aE(4,1)=-0.04884659515311857;
aE(4,2)=-0.17772065232640102;
aE(4,3)=0.8465672474795197;
aE(5,1)=-0.15541685842491548;
aE(5,2)=-0.3567050098221991;
aE(5,3)=1.0587258798684427;
aE(5,4)=0.30339598837867193;
aE(6,1)=0.2014243506726763;
aE(6,2)=0.008742057842904185;
aE(6,3)=0.15993995707168115;
aE(6,4)=0.4038290605220775;
aE(6,5)=0.22606457389066084;
aE(2,1)=2*gama;

aI(3,1)=0.137776;
aI(3,2)= -0.055776;
aI(4,1)=0.14463686602698217;
aI(4,2)=-0.22393190761334475;
aI(4,3)=0.4492950415863626;
aI(5,1)=0.09825878328356477;
aI(5,2)=-0.5915442428196704;
aI(5,3)=0.8101210538282996;
aI(5,4)=0.283164405707806;
aI(2,1)=gama;
for j=2:6
   aI(j,j)=gama; 
end
for j=1:5
   aI(6,j)=b(j); 
end

%This is a file named pro_1221_stable.m.
clear
ARK4_coeffi;
I=eye(6);
B=zeros(6,6);
for j=1:6 
   B(j,:)=b;
end
syms lam_a;
syms lam_d;
R_up=det(I-lam_d.*aI-1i*lam_a.*aE+(lam_d+1i*lam_a).*B);
R_low=det(I-lam_d.*aI-1i*lam_a.*aE);
R=abs(R_up./R_low)-1;
figure(1);
fp1=ezplot(R,[0,15,-80,1]);
view(90,-90);
hold on;
set(gca,'XTick',[0:1:15]);   
set(gca,'yTick',[-80:5:1]); 
axis([0 15 -80 1]);
grid on;
hold on;
title('ARK4');
\end{lstlisting}
	The output is shown below.

	\includegraphics[width=0.45\textwidth]{pro1221_1.png}
	\includegraphics[width=0.45\textwidth]{pro1221_2.png}

	It's clear from the second plot that in the absence of diffusion, the scaled advection eigenvalues should be no more than 4.
	Calculate the possible maxium of $sin\theta_d(1+\frac{2}{3}sin^2\frac{\theta_d}{2})$, $\theta \in (0,\pi)$ by Matlab R2019b.
	\lstset{language=Matlab}
	\begin{lstlisting}
%This is a file named calculate_courant.m.
clear;
syms theta;
y=sin(theta)+2.0/3*sin(theta)*sin(theta/2)*sin(theta/2);
y_max=double(max(subs(y,theta,[0:0.001:pi])))
Cr=4.0/y_max

%This is the output.
>> calculate_courant

y_max =

    1.3722


Cr =

    2.9150
\end{lstlisting}
	Thus, we conclude that stable Courant number for this method should satisfy
	\[ \mu \leq \frac{2.9150}{D}.\]
\end{proof}

\section{13.11 Express $DG$ as linear combination of cell averages to verify $DG\neq L$, and figure out the 1D stencil of $DG$.}
\begin{proof}[Proof]
	In fact, $DG$ can be expressed in
	\begin{align*}
		DG(<\psi>_i) & =\frac{1}{144h^2}(<\psi>_{i+4e^d}-16<\psi>_{i+3e^d}+64<\psi>_{i+2e^d}+16<\psi>_{i+e^d}-130<\psi>_{i} \\
		             & +16<\psi>_{i-e^d}+64<\psi>_{i-2e^d}-16<\psi>_{i-3e^d}+<\psi>_{i-4e^d}).
	\end{align*}
	It's clear that $L\neq DG$ since the stencil width of DG is 9 while that of L is just 5.
	
	For the 1D stencil, consider $E(x)=\int_\xi^{x}\psi$ where $\xi$ is fixed to $x_{i+1}$ to make $E(x_{i+1})=0$. 
	
	Denote $E(x_{i+k}):=E_k$. 
	Then the sentence can be reduced to 
	\begin{align*}
		DG(-E_0)&=\frac{1}{144h^2}(E_5 -17E_4+80E_3-48E_2-146E_1+146E_0+48E_{-1}-80E_{-2}+17E_{-3}-E_{-4}).
	\end{align*} 
	In fact, it's a central difference formula for x=$x_{i+0.5}$.
\end{proof}

\section{13.13 Show $P_E$ is indeed a projection.}
\begin{proof}[Proof]
	\begin{align*}
		P_E^2 & =(I-G(DG)^{-1}D)^2 \\&=I-2G(DG)^{-1}D+G(DG)^{-1}DG(DG)^{-1}D\\&=I-2G(DG)^{-1}D+G(DG)^{-1}D\\&=I-G(DG)^{-1}D=P_E.
	\end{align*}

\end{proof}

\end{document}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 
